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In turbulent thermal convection in cylindrical samples of aspect ratio r = D/L (D is the diameter 
and L the height) the Nusselt number Nu is enhanced when the sample is rotated about its vertical 
axis, because of the formation of Ekman vortices that extract additional fluid out of thermal bound- 
ary layers at the top and bottom. We show from experiments and direct numerical simulations 
that the enhancement occurs only above a bifurcation point at a critical inverse Rossby number 
l/Ro c , with l/Ro c oc 1/r. We present a Ginzburg-Landau like model that explains the existence 
of a bifurcation at finite l/Ro c as a finite-size effect. The model yields the proportionality between 
l/Ro c and 1/r and is consistent with several other measured or computed system properties. 



PACS numbers: 47.27.te,47.32.Ef,47.20.Bp,47.27.ek 

Turbulence, by virtue of its vigorous fluctuations, is 
expected to sample all of phase space over wide parame- 
ter ranges. This viewpoint implies that there should not 
be any bifurcations between different turbulent states. 
Contrary to this, several cases of discontinuous transi- 
tions have been observed recently in turbulent systems 
[T]-[3]. When they occur, they are likely to be provoked 
either by changes in boundary conditions or boundary- 
layer structures, or by discontinuous changes in the large- 
scale structures, as a parameter is varied. 

Recently some of us [11 [5] reported on the effect of 
rotation about a vertical axis at a rate Q on turbulent 
convection in a fluid heated from below and cooled from 
above (known as Rayleigh-Benard convection or RBC; 
for recent reviews, see [6HH])- For a cylindrical sample 
of aspect ratio F = D / L = 1.00 (D is the diameter and 
L the height) a supercritical bifurcation was found, both 
from experiments and from direct numerical simulation 
(DNS) of the Boussinesq equations of motion. At a finite 

as expressed by the inverse Rossby number 1/Ro ex Q 
(to be defined explicitly below), there was a sharp tran- 
sition from a state of nearly rotation-independent heat 
transport (as expressed by the Nusselt number Nu to 
be defined explicitly below) to one in which Nu was en- 
hanced by an amount 5Nu(l/Ro). This is illustrated by 
the data shown in Fig. [I] The increase of Nu was at- 
tributed to Ekman pumping [4) \5\ 19HT5], i.e. to the for- 
mation of (cyclonic) vertical vortex tubes ("Ekman vor- 
tices"), which extract and vertically transport additional 
fluid from the boundary layers (BLs) and thereby en- 
hance the heat transport. The bifurcation was located 
at a critical value l/Ro c ~ 0.40 [5]. The reason for the 
existence of the bifurcation at l/Ro c > hitherto had 



not been understood. While such bifurcations are com- 
mon near the onset of RBC in the domain of pattern 
formation [16], their existence in the turbulent regime 
implies a paradigm shift. 

In this Letter we report on further experiments for 
samples with F = 2.00, 1.00, and 0.50 which (i) all show 
bifurcations between different turbulent states and (ii) 
reveal that l/Ro c varies approximately in proportion to 
1/T. We offer an explanation of these and other phenom- 
ena in terms of a phenomenological Ginzburg-Landau like 
description which predicts a finite-size effect upon the 
vortex density A. 

We assumed that the relative Nusselt enhancement 
£Nu(l/Ro)/Nu(0) is proportional to the average A of A 
over a horizontal cross section of the sample near the 
BLs. Consistent with the DNS that we report here, 
we assumed that A vanishes at the sample side wall. 
For the infinite system the model predicts that A, and 
thus 5Nu/Nu(0), increases linearly from zero starting at 
1/Ro = 0. For the finite system the model gives a thresh- 
old shift proportional to 1/T as found in the experiment 
and by DNS. The shift is predicted to be followed by a lin- 
ear increase of A in proportion to (1/Ro) — (l/Ro c ) which 
yields JNu/Nu(0) = Si (r) (1/Ro - l/Ro c ). The model 
gives an initial slope Si(F) that decrease with decreasing 
T, again consistent with DNS and measurements. From 
DNS we show that A decreases to zero near the side wall 
over a length that is consistent with an estimate of a 
healing-length £ based on the model. Thus, we found 
consistency between the model predictions and all prop- 
erties that we were able to either measure or compute 
from DNS. 

Before proceeding, we define the relevant dimension- 
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FIG. 1. (Color online) The Nusselt number Nu(l/Ro), nor- 
malized by Nu(0) without rotation, as a function of the inverse 
Rossby number 1/Ro. Data are from experiments unless men- 
tioned otherwise, (a): r = 0.50 and Pr = 4.38; Rayleigh num- 
bers Ra = 9.0 x 10 9 (solid circles, red online), 1.8 x 10 10 (open 
circles, purple online), and 3.6 x 10 10 (solid squares, blue on- 
line), (b): T = 1.00. Main figure: Pr = 4.38, Ra = 2.25 x 10 9 
(solid circles, red online), 8.97 x 10 9 (open circles, purple 
online), and 1.79 x 10 10 (solid squares, blue online). Inset: 
Pr = 6.26 and Ra = 2.73 x 10 8 . DNS: open squares (black 
online). Experiment: solid circles (red online), (c): r = 2.00 
and Pr = 4.38; the data are for Ra = 2.91 x 10 8 (open cir- 
cles, red online, DNS; solid circles, purple online, experiment), 
5.80 x 10 8 (open squares, blue online), and 1.16 x 10 9 (solid 
squares, black online). Note the different vertical scale for (c) 
compared to (a) and (b) . The small vertical lines indicate the 
locations of the bifurcation points. 



less parameters. The inverse Rossby number is given 
by 1/Ro = (2Sl)/(y/PgAT/L), were ft is the rotation 
rate in rad/s, (3 the isobaric thermal expansion coeffi- 
cient, AT the temperature difference between the bot- 
tom and top plate, and g the gravitational acceleration. 
The Rayleigh number is Ra = (PgATL 3 )/(vK), were 
k, and v are the thermal diffusivity and the kinematic 
viscosity respectively. The Nusselt number is given by 




FIG. 2. The critical inverse Rossby number l/Ro c as a func- 
tion of the inverse aspect ratio 1/T for Pr = 4.38 and different 
Ra (see Fig.[T]). The line is a fit of a quadratic equation with 
the constant coefficient set to zero. 



Nu = (QL)/(ATX) where Q is the heat-current density 
and A is the thermal conductivity. Finally, the Prandtl 
number is Pr = v j k. 

In Fig. [I] we show experimental and numerical data [? 
] for Nu(l/Ro)/Nu(0) as a function of 1/Ro for several 
values of Ra. From top to bottom the three panels are 
for T = 0.50, 1.00, and 2.00, respectively [? ]. One sees 
that there is considerable structure even below the bifur- 
cation, particularly at the larger Ra. To our knowledge 
the origin of this structure is not known in detail. One 
sees that there are clear breaks in the curves, e.g. for 
T = 1.00 (Fig. [TJ)) at 1/Ro ~ 0.4, indicating the bifurca- 
tion to a different state. The location of this transition 
is within our resolution independent of Ra. 

In Fig. [2] we plotted all available data for 1 /Ro c for 
Pr = 4.38 (and different Ra) as a function of 1/T. The 
line shown there is a fit of 
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to the data. Its coefficients are a = 0.381 and b = 0.061. 
One sees that the data are consistent with an initial lin- 
ear increase from zero of l/Ro c with 1/T, with a small 
quadratic contribution becoming noticeable as 1/T be- 
comes larger. 

In order to understand the T dependence of l/Ro c , 
we studied the vortex statistics using data obtained from 
DNS. We used the so-called Q-criterion |TTJ [T7HT9] to 
determine the percentage A of the horizontal area that 
was covered by vortices. Using this criterion, implies 
that the quantity Q2D [3 HO], which is a quadratic form 
of various velocity gradients, was calculated in a plane 
of fixed height. An area is then identified as "vortex" 
when Q 2D < -0.1 (\Q 2 d\) 

v , were (\Q2d\)v is the volume- 
averaged value of the absolute values of Q2D [5 - The 
result of this procedure is shown for different 1/Ro in 
Fig. [3] for Pr = 6.26. In Fig. [4] we plot I as a function 
of 1/Ro at the edge of the kinetic BL (which depends 
on Ro, see [21]) and at the fixed distance 0.023L (the 
kinetic BL thickness without rotation) from the plates. 




FIG. 3. The vortices as identified by the Q^d criterion [51 [20] 
for Ra = 2.73 x 10 8 , Pr = 6.26, and r = 1. a) 1/Ro = 2/3, b) 
1/Ro = 1, c) 1/Ro = 1.54, and d) 1/Ro = 3.33. The vortex 
area increases with increasing 1/Ro. This trend is quantified 
in Fig. [4] 
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FIG. 4. Fraction A of a horizontal slice covered with vortices 
at the edge of the kinetic BL (circles) and at a distance 0.023L 
(the kinetic BL thickness without rotation, squares) from the 
plates as a function of 1/Ro for Ra = 2.73 x 10 8 , Pr = 6.26, 
and r = 1. The vertical dashed line indicates the bifurca- 
tion point at l/Ro c and the horizontal dash-dotted line is a 
background vorticity level Aq present even below l/Ro c . 



Although there is quite a bit of scatter, the data are con- 
sistent with a linear increase of A for 1/Ro > l/Ro c , with 
a small constant background A = Aq below l/Ro c . The 
azimuthally averaged vortex density (A)^ is given in Fig. 
[5] It shows that the Ekman vortices are inhomogeneously 
distributed: While in the bulk their fraction is roughly 
constant - the variations are (presumably) due to insuf- 
ficient statistics - there are almost no vortices at all close 
to the side wall, signaling a strong boundary effect. 

In an effort to understand the existence of a finite on- 
set (see Fig. [I]) of the Ekman- vortex formation and the 
dependence of the critical inverse Rossby number on T 




FIG. 5. Azimuthal average (A)<f,(r/L) of the vortex density 
A for Ra = 2.73 x 10 8 , Pr = 6.26, r = 1, and 2.22 < 1/Ro < 
3.33. In total, the statistics is based on 8 snapshots. The 
dashed line is the uniform case. The density approaches zero 
close to the side wall. 



(see Fig. to elucidate the linear rise and initial slope 
of Nu(l/Ro) above onset (see Fig. [I]), and to explain the 
rapid decrease of (A)(p(r/L) near the wall (see Fig.[5|, we 
propose a phenomenological Ginzburg-Landau like model 
for the local vortex density 



A = (1/Ro 2 ) A - gA 3 + £ 2 V 2 A 



(2) 



Here A is the time derivative of A. We chose the coef- 
ficient of the linear term as 1/Ro 2 because for the time 
independent infinitely extended spatially uniform system 
it yields a stable solution A = g -1 / 2 (l/Ro) which implies 
a vortex density proportional to the rotation rate. The 
term with V 2 A represents the lowest-order term of a gra- 
dient expansion since terms proportional to VA would 
lead to an unphysical propagating mode. 

When spatial variations are allowed, the ground state 
A = can be shown to be stable (i.e. to have a growth 
rate a < 0) to disturbances with wave vector k when 
1/Ro falls below a neutral curve given by 



l/Rx>o(fc) =£ok . 



(3) 



For the finite system it is necessary to introduce ap- 
propriate boundary conditions. Here we shall consider 
a one-dimensional system over the range — T/2 < x < 
T/2 for simplicity and illustrative purposes. The two- 
dimensional system with circular boundaries and no az- 
imuthal variation was treated in detail in Ref. [22] and 
yields the same result for l/Ro c . Since there can be no 
vortices at the side wall of the sample (see Fig. [5] where 
we verified this based on the numerical data), we chose 
A(-T/2) = A(T/2) = 0. For the wave number fc of the 
lowest mode this yields ko = n/T. This in turn gives 

JL = 1 = <o 

Ro. 



(4) 



Roo(^o) r ' 

Thus, consistent with the data in Fig.|2j the model yields 
the proportionality between l/Ro c and 1/r. We note 
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that the curvature indicated by the quadratic contribu- 
tion to the fit Eq. [I] can be accommodated easily by 
higher-order gradient terms in Eq. [2] Comparison with 
experiment (see Eq. [I]) gives £o = cl/tt — 0.121. 

To elucidate the rapid decrease of (Aj^r/L) in Fig. [5] 
near r/L = 0.5, we consider Eq. [5] for a semi-infinite 
system over the range — oo < x < 0.5 with the boundary 
condition A(x = 0.5) = 0. It yields the solution 

A(x) = (Ro 2 g)~ 1/2 tanh ((0.5 - x)/£) (5) 

with 

£ = v^CoRo . (6) 

Thus, near the boundaries, the model predicts that the 
amplitude A(x) of the one-dimensional model, and thus 
to a good approximation also the azimuthal average 
(A) ( f ) (r/L) in Fig. [5| should "heal" to its bulk value over 
a length £. Using a representative Ro ~ 0.4 for Fig.|5j we 
estimate from Eq.|6]that £ ~ 0.07. This is roughly consis- 
tent with the rapid variation of (A)^(r / ' L) near r/L = 0.5 
seen in Fig. [5] 

Above the bifurcation the model yields [22] 



Thus 




which is consistent with the data in Fig. [T] Numerical so- 
lutions of the amplitude equation of Ref. [22] have shown 
that the re- normalized coefficient g in Eq.[7]is larger than 
g in Eq. [2] Thus, the initial slope of A above l/Ro c and 
Si (r) in Eq. [8] are reduced by the finite size of the sys- 
tem. The decrease of Si(r) with decreasing T that can 
be seen in Fig. [I] is also consistent with the model. 

At constant T, Si depends slightly on the Prandtl num- 
ber. This suggests that the nonlinear coefficient g in 
Eq. [2] and thus g in Eq. [7| is dependent on Pr. Simi- 
larly, the bifurcation point l/Ro c depends slightly on Pr. 
This is accommodated in the model Eq. [2] by a slightly 
Prandtl-dependent length scale £o- 

It remains to be seen whether the phenomena re- 
ported and explained here in terms of a finite-size effect 
have analogies in bifurcations between turbulent states 
in other systems [TH2]. The more general lesson which 
is learned is that the Ginzburg-Landau approach, which 
has been so versatile to understand the spatio-temporal 



dynamics of patterns, can also be useful in understanding 
the remarkable bifurcations between turbulent states. 
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